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ABSTRACT 

We model the formation and evolution of galaxy clusters in the framework of an 
extended dark matter halo merger-tree algorithm that includes baryons and incorpo- 
rates basic physical considerations. Our modified treatment is employed to calculate 
the probability density functions of the halo concentration parameter, intracluster 
gas temperature, and the integrated Comptonization parameter for different cluster 
masses and observation redshifts. Scaling relations between cluster mass and these 
observables are deduced that are somewhat different than previous results. Modeling 
uncertainties in the predicted probability density functions are estimated. Our treat- 
ment and the insight gained from the results presented in this paper can simplify 
the comparison of theoretical predictions with results from ongoing and future cluster 
surveys. 
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1 INTRODUCTION 

Formation of galaxy clusters is of central importance for un- 
derstanding the evolution of the large scale structure (LSS) 
of the universe. Statistical properties of clusters - deduced 
from cluster optical, X-ray, and SZ surveys - can be used to 
determine the basic cosmological parameters - such as the 
matter density, the normalization (and spectrum) of primor- 
dial density fluctuations, and the dark energy equation of 
state - independently of other methods (CMB power spec- 
trum, galaxy surveys, etc). This can only be achieved if the- 
oretical tools are developed for a quantitative description of 
the non-linear hierarchical growth of clusters. Hydrodynam- 
ical numerical simulations of clusters are currently the most 
versatile tool for establishing statistical properties of clus- 
ters, which are commonly specified in terms of mass func- 
tions and scaling relations between their intrinsic properties. 

In order to use cluster surveys as cosmological probes it 
is essential to know how their intrinsic properties - such as 
formation time, mass and redshift - affect their integrated 
statistical properties. Clusters are commonly described as 
virialized spherical systems whose masses are dominated by 
dark matter (DM) and with isothermal intracluster (IC) gas 
as their main baryonic mass component. Cluster statistical 
relations are usually expressed in terms of the mass and the 
redshift of observation; for example, the resulting scaling of 
the gas temperature is 
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where E(z) = H (z) / Ho, the Hubble parameter in units of its 
present value and Ay(z) is the overdensity at virialization. 

However, this simple description may not be sufficiently 
adequate for the description of real clusters. Observations 
and numerical simulations indicate that DM density profiles 
are shallower than isothermal at small radii, and steeper 
than isothermal at large radii, and are characterized by a 
scaling radius r s which marks the transition between these 
two regions (Navarro, Frenk & White 1995, 1996). Although 
it seems reasonable to assume that the scaling relations are 
not much affected by the details of the cluster structure, the 
standard description is not sufficiently accurate since the 
model parameters - such as r s - depend on the cluster mass 
and redshift. Indeed, numerical simulations show that the 
formation history of clusters affects the DM scaling radius, 
such that clusters are systematically denser the earlier they 
formed. Cluster mass concentration is usually quantified in 
terms of the parameter c, defined as the ratio of the cluster 
virial radius to the scaling radius, c = R v /r s . Thus, clus- 
ters that form earlier have larger concentration parameters 
(Wechsler et al. 2002); more generally, other physical prop- 
erties are also influenced by the cluster formation history. 

The assumption of hydrostatic equilibrium is obviously 
an approximation, since clusters, the largest bound systems, 
are still forming through mergers of subclumps and accre- 
tion. Mergers disrupt the state of thermal equilibrium; dur- 
ing some merger events IC gas temperature and X-ray lumi- 
nosity are boosted up by factors of up to 3 and 10, respec- 
tively, as was demonstrated in a series of numerical simu- 
lations by Ricker and Sarazin (2001). More generally, when 
deriving the standard scaling relations one usually ignores 
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the formation history of the cluster. In the ACDM model 
the LSS formed hierarchically, through consecutive merg- 
ers of smaller structures, and thus clusters with roughly the 
same mass and redshift may have very different merger and 
formation histories. Indeed, both observations and simula- 
tions reveal scatter in the mass-observable relations at a 
level of ~ 10% — 20%, which is partly due to the differ- 
ent formation histories (e.g., Wechsler et al. 2002, Vikhlinin 
et al. 2009a). In order to use these relations to determine 
the cluster mass function from surveys it is important not 
only to understand the exact scaling relations, but also to 
quantify the amount of scatter. Furthermore, uncertainties 
in the mass-observable relation reduce the precision with 
which cosmological parameters can be determined (Lima & 
Hu 2005; Cunha & Evrard 2009). 

In current analyses of cluster surveys the scatter in 
the mass-observable relations has been partly accounted for. 
Some of this scatter is clearly due to the different dynamical 
state of the clusters in the survey; unrelaxed systems are 
sometimes excluded from the analysis. Several studies have 
indicated that the scatter in the scaling relations is reduced 
if IC gas temperature is measured by excluding the cluster 
central region (which is significantly affected by radiative 
cooling), but the intrinsic scatter in the scaling relations is 
more difficult to estimate. For example, motivated by insight 
from simulations, Vikhlinin et al. (2009b) used a constant 
value of 20% as an estimate for the Tx — M scatter. How- 
ever, the amount of scatter may depend on the mass and 
redshift of the cluster. 

A meaningful comparison with observational data usu- 
ally necessitates knowledge of the full probability distribu- 
tion function (PDF) rather than just the scaling with mass. 
An example is the case of very large values of the concen- 
tration parameter and Einstein radius for several clusters 
(Broadhurst et al. 2008, Zitrin et al. 2009). Based on results 
from N-body simulations, the probability of observing clus- 
ters with the very high measured values was found to be 
very low, amounting to a 4 — a discrepancy with the ACDM 
predictions. The concentration parameter PDF is a crucial 
component in this analysis (Sadeh & Rephaeli 2008). It is 
clear that a better understanding of the origin of this PDF 
is required. 

The impact of the formation history on the cluster prop- 
erties was previously investigated using a series of hydrody- 
namical simulations (Ricker and Sarazin 2001, Randall and 
Sarazin 2002, Wik et al. 2008). These authors used simu- 
lations of pairs of merging clusters and analysed the im- 
pact of recent mergers on IC gas temperature, luminosity, 
and Comptonization parameter. They found that all these 
quantities are boosted for a relatively long time following 
a merger, and calculated the effect of this boost on cos- 
mological parameter estimation. These analyses did not in- 
clude the full formation history, only the latest merger, and 
relied on a small number of simulated clusters with differ- 
ent masses. Voit and Donahue (1998) showed that the tem- 
perature evolves less rapidly with mass than in the stan- 
dard analysis when the recent formation approximation is 
relaxed. They assumed gradual mass accretion throughout 
the cluster history, and a one-to-one correspondence between 
the temperature and the virial energy of the cluster. 

The statistics of DM halo concentrations and their de- 
pendence on the halo formation history were investigated 



using N-body simulations (Bullock et al. 2001, Wechsler et 
al. 2002, Neto et al. 2007, Gao et al. 2008, Duffy et al. 2008). 
This approach provides PDFs of the concentration parame- 
ter which account for different formation histories of differ- 
ent halos. Results of these works can then be used to infer 
the intrinsic scatter in other cluster observables. However, a 
theoretical approach that is complementary to N-body sim- 
ulations is needed in order to fully understand the impact of 
the cluster formation history on its properties. The reliabil- 
ity of the statistical analysis of N-body simulations depends 
on the simulation volume, which is limited by computational 
constraints. This can be a severe problem if one is interested 
in high-mass clusters, which are relatively rare systems. For 
example, the Millennium Simulation (MS, Springel et al. 
2005), the largest cosmological N-body simulation to date, 
which follows N = 2160 3 particles in a periodic box of 
L = 500/i." 1 Mpc on a side, contains less than 800 relaxed 
halos with masses above M200 = 1.3 • 1O 14 /i _1 M0, and just 
8 relaxed halos with masses above M200 = 7.5 ■ lO 14 ft _1 M0 
(Neto et al. 2007, hereafter N07). Moreover, comparison 
between different simulations is difficult because each uti- 
lizes different cosmological parameters and halo finding al- 
gorithms. This difficulty is illustrated by the fact that some- 
what different mass functions are predicted by different N- 
body simulations (Jenkins et al. 2001, Sheth & Tormen 2002, 
Tinker et al. 2008), most likely reflecting the different for- 
mation histories predicted by these simulations. 

Predicting the full PDFs of the relevant cluster param- 
eters in the context of a theoretical model that can be read- 
ily implemented, would allow a more meaningful statistical 
analysis and the ability to quantify the impact of uncertain- 
ties in the values of cosmological and cluster parameters on 
the main observables deduced from large scale surveys. In 
this paper we develop a model of cluster formation using an- 
alytically computed DM merger trees, with which we trace 
the formation of clusters through major episodal mergers 
and continuous accretion. We show that our approach pro- 
vides an improved physical description in comparison with 
what can be obtained from standard scaling relations. 

This paper is organized as follows. In Section [2] we de- 
scribe our model of cluster formation. The results of our 
generalized merger-tree treatment are presented in Section 
[3] and further discussed in Section |4] Throughout the paper 
we use the following cosmological parameters: Q m = 0.25, 
Q A = 0.75, H = 73 km/s/Mpc, a 8 = 0.8. 



2 METHODOLOGY 

Our description of the growth of galaxy clusters is based 
on merger trees of DM halos as described numerically by 
the modified GALFORM code (Cole et al. 2000; Parkin- 
son, Cole & Helly 2008), which was successfully employed 
to construct semi-analytic models of galaxy formation. The 
algorithm implements the excursion set formalism, a key as- 
pect of which is the conditional mass function: the fraction 
of mass /(M1IM2) from halos of mass M2 at redshift 22 that 
consisted of smaller halos of mass Mi at an earlier redshift 
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where <7; = a 2 (Mi) is the variance of the linear perturbation 
field smoothed on scale Mi, and Si is the critical density for 
spherical collapse at redshift z;. The original GALFORM al- 
gorithm is consistent with the Press- Schechter mass function 
(Press & Schechter, 1974; PS), in the sense that if a grid of 
trees is rooted at z = 0, weighted by their mass abundance 
according to the PS mass function, then the mass function 
at higher redshifts again corresponds to PS mass function. 
The conditional mass function is used to calculate the mean 
number of progenitors of mass Mi at redshift zi + dzi of a 
halo of mass M2 at 22 = zi : 



dN _ J_ / df 

dMi ~ Mi I dzi 



M2 
Mi 



dzi 
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The modified GALFORM algorithm was obtained by mak- 
ing the substitution 

dN dN 



dMi dMi 



G(ffl/(T2, $2/(72) 



(4) 



where G is referred to as a perturbing function, to be cal- 
ibrated by comparison with N-body simulations. Parkin- 
son et al. (2008) showed that by fitting the outcome of 
the algorithm to the results of the MS they obtained halo 
abundances which are consistent with the Sheth-Tormen 
mass function (Sheth & Tormen 2002). Thus, the perturb- 
ing function G expresses the uncertainty in the choice of 
the correct mass function. In this work we use the fol- 
lowing parametrization of the perturbing function: G = 
Go (o~i/a 2 )' yi (#2/o"2) 72 , with the parameters Go = 0.57, 
71 = 0.38, 72 = -0.01 taken from Parkinson et al. (2008). 

Starting with the specified mass and redshift, the algo- 
rithm proceeds back in time, checking after each timestep 
whether there was a merger; if so, the masses of the merged 
halos are drawn from the distribution ©. Halos with masses 
below some resolution limit M res are not resolved, and are 
accounted for as continuously accreted mass. Further details 
on the GALFORM algorithm can be found in Parkinson et 
al. (2008). 

In the following sections we describe our model of clus- 
ter formation. Since the DM is the dominant mass com- 
ponent of clusters we first study the formation history of 
cluster-sized DM halos, and then add the IC gas compo- 
nent, study its properties and related scaling relations. 



2.1 Modeling the formation of dark matter halos 

The calculation begins with the construction of a merger tree 
for a given final halo mass and redshift. For each tree, we 
use only the major merger events, that is only those merg- 
ers for which M> /M < < q, whose value is to be determined. 
The rationale behind this choice is the assumption that the 
properties of a cluster are largely determined by the violent 
merger events, during which DM and gas settle in the modi- 
fied potential well when equilibrium is reestablished. On the 
other hand, slow accretion of material on the outskirts of the 



cluster, as well as mergers with low-mass systems (which are 
treated in an identical manner in this work) do not cause re- 
distribution of the cluster components and do not strongly 
affect the physical conditions in the cluster center. Thus, we 
follow the major merger events in the cluster history, and 
refer to all other processes of mass growth as slow accre- 
tion. A typical value of the major merger ratio is q = 10; 
other values will also be considered. The trees were gener- 
ated up to the redshift z — Zf in which depends on the mass 
resolution of the merger tree algorithm, M res , as discussed 
below. 

The next step is to calculate the density profile of each 
halo in the tree. To this end, the tree is traced down, begin- 
ning with the smaller masses; for each merger event we use 
energy conservation to calculate the density profile of the 
merged DM halo. For simplicity we assume NFW (Navarro, 
Frenk & White 1995) profiles for all the halos at all times: 



Pdm(t) 



x(l + x) 



(5) 



where x = r/r s = rc/R v is the radial distance expressed 
in terms of NFW scale radius r s , c is the concentration pa- 
rameter, and p s is the mass density normalization constant. 
Each halo is completely characterized by its mass, redshift 
and concentration parameter, where halo mass is defined 
as the mass within the virial radius, M = ^-R^pcAy- All 
masses are assumed to be in equilibrium, which is presumed 
to be attained relatively quickly after each merger event, but 
we do not assume the halos are completely virialized within 
the virial radius. In fact, the virial ratio 2T/\W\ approaches 
1 at the virial radius only for very large c, while for com- 
monly deduced values of c this ratio is slightly larger than 
1 at the virial radius (Cole & Lacey 1996; Lokas & Mamon 
2001). 

The relation between mass and virial radius depends on 
redshift, which we identify as the halo formation redshift, 
Zf, defined as follows. If the halo has undergone a major 
merger event, we take Zf to be the redshift of the last major 
merger. This choice is consistent with the main assumption 
that major mergers largely determine the physics of the halo. 
For halos that did not experience major mergers at all, and 
were, according to our interpretation, entirely assembled by 
minor mergers and continuous accretion, we take Zf to be 
the redshift at which half of the halo mass had assembled. 
Difficulty in determining zj by this prescription is encoun- 
tered only when the branch of the merger tree terminates 
when the halo still has more than half of its mass (recall 
that the tree is evolved backwards in time). In GALFORM, 
the branch is terminated in two cases: either Zfi„ is reached, 
or the mass of the halo falls below the resolution mass M res . 
Thus, in order to describe the formation of the smallest ha- 
los that constitute the tree, Zfi„ should clearly be chosen 
well above the expected formation redshift of halos of mass 
Mres- By making this choice we ensure that almost all of the 
halos in the tree assemble half of their mass before z/; n is 
reached (that is, later in time), and their formation redshift 
can be traced back by the tree. 

To start from the earliest halos in the tree and move 
forward in time requires specifying their concentration pa- 
rameters. These are adapted from a fit to a set of N-body 
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simulations by Bullock et al. (2001) 
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where zj and z (, s are the redshifts of formation and obser- 
vation, respectively. This choice is motivated by the finding 
of Wechsler et al. (2002) that the concentration parameter 
scales as c ~ (1 + z/)/(l + z ts), although note that their 
definition of the formation redshift is slightly different than 
the one adopted here. Although this choice for the initial 
c(M, z) is somewhat arbitrary, its particular form does not 
significantly influence the results. 

For each merger event, we calculate the total energy 
of the system before merging, which depends on the con- 
centration parameters of the merging halos, and deduce the 
concentration parameter of the merged halo from simple en- 
ergy conservation arguments, motivated by a cluster merger 
model by Sarazin (2002). The total energy of a system of 
two halos before merging is: 



E tot ,xi = E{M X ) + E{M 2 ) + U 12 



(7) 



where E(Mi) are the total energies (potential and kinetic) 
of each halo and U12 is the gravitational energy of the two 
halos at the point of their largest separation, when they have 
just become bound and their relative velocity was negligible: 



Ui 



GM1M2 
d 



(8) 



The distance do, at which the halos became bound, is 
roughly the mean distance between halos with masses Mi 
and M2 that reside in an overdense region with a scale that 
corresponds to the final mass Mf = Mi + M2 . We take this 
distance to be do = red, where d = Ri + R2, and adopt re = 5 
as a fiducial value for all halos. This corresponds to a typ- 
ical distance of several Mpc and a typical relative velocity 
of several hundreds to a few thousands of km/s, depending 
on the masses of the merging clusters, which is in accor- 
dance with the initial conditions of hydrodynamical cluster 
merger simulations (Ricker & Sarazin 2001, McCarthy et al. 
2007, Lee and Komatsu 2010). Very high values of re produce 
unrealistically large initial separations and large relative ve- 
locities, while very low values of re lead to very small initial 
distances, small relative velocities, and consequently, greatly 
reduced total energies, which eventually result in very high 
concentration parameters of the final halo. In other words, 
re was chosen so as to yield realistic values of both the ini- 
tial separation and the final concentration parameter. The 
dependence of the results on re is discussed below. 

After the two halos merge, the resultant halo accretes 
matter, so when in turn it merges to form a larger halo, it 
has more mass than the sum of the masses of its progenitors. 
We account for the energy of the accreted matter in a very 
approximate way as follows. Given the masses of the two 
progenitor halos, Mi and M2, and the final mass of the halo, 
M p , the total accreted mass is AM = M p - (Mi + M 2 ) = 
M p — Mf. Numerical simulations of galaxy-sized halos (e. g. 
Wang et al. 2010) seem to indicate that the accreted material 
is distributed in the halo outer region. This is quite likely 
the case also in cluster-sized halos, so we can estimate the 



energy due to the accreted mass by writing 



Uacc — 



G(Mi + M 2 )AM 



R 



f 



(9) 



where Rf is the virial radius of the halo with mass M/ = 
Mi + M2 that formed just after the merger. In evaluating 
Uacc we assume that the accreted mass constitutes a rela- 
tively small fraction of the final halo, and that equation ([9} 
provides a simplified description of the accretion process. 

We then have for the total energy of the system prior 
to merging 



Etotal — EtotA2 + Uacc 



(10) 



This energy is attributed to the resulting merged halo; we 
assume no mass is lost in the process. By equating E to tai 
with the gravitational and potential energy of the resulting 
halo, which depend on its concentration parameter, we can 
deduce the latter. This process is repeated for each halo in 
the tree, until arriving at the bottom - the most massive 
halo. At the end of this process we obtain the concentration 
parameter for the given mass and for one tree, i.e. one re- 
alization of the halo history. Generating a large number of 
trees gives an estimate of the PDF of c(M). 



2.2 Modeling the intracluster gas 

In our modeling of IC gas we assume that it constitutes a 
small fraction of the total cluster mass, and that it does not 
significantly affect the evolution of the cluster. We assume 
that the gas has a polytropic equation of state with an adi- 
abatic index V, such that the gas density and pressure are 
related through: 



Po(p/f>o) 



(11) 



The solution of the equation of hydrostatic equilibrium 
for a polytropic gas inside a potential well of a DM halo with 
an NFW profile is (Ostriker, Bode, & Babul 2005): 



p{x) = p 



1 



B 



1 + n 



ln(l + x) 



where n = (T — 1) 1 and B is given by: 



B = 



4,-KGp s r 2 pm p 
fcsTo 



(12) 



(13) 



and fim p is the mean molecular weight. Thus, B depends on 
the concentration parameter through p s and r s (see equation 
(0). The temperature profile is then 



T(x) = T 



1 - 



B 



1 + n 



1 - 



ln(l- 



(14) 



As a boundary condition we assume that the gas pres- 
sure at the virial radius obeys P ga s = fgPDM where f g is the 
gas mass fraction and Pdm = Pdm& 2 (Ostriker et al. 2005). 
We obtain a 2 , the DM (3D) velocity dispersion, by solving 
the Jeans equation for the NFW potential. We expect this 
particular choice for the boundary condition to have a minor 
influence on the results, as discussed in Section [4] 

For each merger tree we obtain the concentration pa- 
rameter and the virial radius of the final halo, as described 
above. Taking the observationally deduced value for the adi- 
abatic index, F = 1.2, and imposing the above boundary 
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condition to obtain the constant B fully determines the tem- 
perature profile. By assuming a specific gas mass fraction we 
also obtain the full density profile. Repeating this procedure 
for a large number of trees provides an estimate of the PDFs 
of the various physical parameters as a function of the clus- 
ter mass and the redshift of observation. 

In order to derive scaling relations we define the mean 
cluster temperature as the emission-weighted value 



jA(T)p 2 g TdV 
fA(T)p>dV ' 



(15) 



where integration is over the cluster volume, and the temper- 
ature dependence of the cooling function is approximately 
A(T) oc VT. Similarly, X -ray luminosity is approximated 
by: 



Lx = 



( Pa 



\/im p 



A(T)dV. 



(16) 



Knowledge of the gas temperature allows also calcula- 
tion of the Comptonization parameter, y, defined as 



V = 



( k B T e 
\m e c 2 



n e OTdl, 



(17) 



where the integral is taken along the line of sight to the 
cluster, <tt is the Thomson scattering cross section (and the 
electron temperature T e is assumed to equal the gas tem- 
perature). The measured SZ intensity (change) is approx- 
imately proportional to the integrated Y-parameter, given 
by the integral of y over the angle the cluster subtends on 
the sky 



Y 



ydO, 



(18) 



3 RESULTS 

3.1 Model parameters 

The model contains several physical parameters, and two 
numerical (code-specific) parameters - the number of tree 
realizations used to estimate the PDFs and the resolution 
mass of the merger tree. We shall discuss the latter parame- 
ters here and defer the discussion of the physical parameters 
to subsequent sections. 

A key objective of our model is to determine the PDF 
of the concentration parameter and IC gas temperature by 
generating a large number of merger tree realizations N. 
We have found that taking N = 5 ■ 10 4 is sufficient to obtain 
convergent results - taking larger N does not change the 
PDF by more than a fraction of a percent. In what follows, 
we show histograms of 50 binned values obtained from 5 TO 4 
merger trees. 

As noted earlier, a resolution mass needs to be selected 
for each tree. This mass is the smallest building block used in 
the tree. By sampling different values of the resolution mass, 
we find that M rea has to be at least 3 orders of magnitude 
below M, the final mass for which the tree is built, while 
taking smaller values of M res does not affect the results: for 
the mass range we consider, the mean value of c changes 
by no more than 1% when the resolution mass is lowered 
from M res = 10~ 3 A'/ to M res = 10~ 4 Af. The value of M res 
determines z/i n , as discussed above. 
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Figure 1. Concentration parameter probability distribution 
function (PDF) as predicted by the merger-tree model (bars), the 
lognormal fit to this PDF (solid line), a distribution extracted 
from the MS by N07 (dot-dashed line) and the prediction of 
SR08 (dashed line). Upper panel: M = 4 X 10 14 /i _1 M Q for the 
merger tree model and a corresponding range of M200 = 10 14,25 — 
10 i4.75 fc -l MQ from the Mg Lower pane i M = 6 X 10 13 h~ 1 MQ 

for the merger tree model and a corresponding range of M200 = 
10 13.63 _ iQiS'SSft-iMQ from the MS. 



3.2 PDF of halo concentration 

The basic outcome of the model is the concentration pa- 
rameter of the DM halo at a given redshift of observation. 
Each merger tree results in a slightly different concentra- 
tion parameter, which depends on the particular structure 
of the merger tree. Thus, in the limit of a large number of 
tree realizations the distribution of formation histories pro- 
vides a PDF of the concentration parameter. The PDF of 
c for M = 4 x lO 14 ft _1 M at z = is shown in Figure 
[1] (upper panel). A log- normal distribution provides a rea- 
sonable fit, with < log 10 c >= 0.706 and o"i oglo c = 0.106. 
The width of the distribution is comparable with the value 
obtained by N07 for a population of relaxed halos in the 
corresponding mass range M 20 o = 10 14 25 - lO 14 75 ft _1 M 
seen in the MS: < log 10 c >= 0.663 and 0"i oglo c = 0.092. 
The distribution for a lower mass of M = 6 x lO 13 fc _1 M 
is shown in the lower panel, along with a corresponding 
distribution for halos in the MS in the mass range of 
M200 = 10 1363 - 10 1388 /i _1 Mq. For this mass we obtain 
< log 10 c >= 0.758 and <rio gl0 c = 0.106 from the log-normal 
fit, compared with < log 10 c >= 0.744 and o"io gl0 c = 0.094 
for the halos in the MS. We note that the mass correspon- 
dence is only approximate, since the relation between M200 
and M v for a given halo depends on its concentration pa- 
rameter. 
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Figure 2. PDF of the concentration parameter from the merger-tree model for M = 10 15 h~ 1 Mq at z = with different model 
parameters, as discussed in the text. Shown in the Left upper panel is the dependence on the mass function through the perturbing 
function G (see equation JD), the dependence on the initial conditions for e(M, z) for the earliest halos in the tree in the Right upper 
panel, the dependence on the major merger parameter q = M>/M< in the Left lower panel, and the dependence on re = do/d (which 
parametrizes the initial distance between halos that are about to merge) in the Right lower panel. 



For comparison, a semi-analytical calculation adapted 
from Sadeh and Rephaeli (2008; SR08) is also shown. This 
latter treatment was based on an analytical distribution of 
formation times and a relation between the formation time 
and concentration parameter deduced from numerical simu- 
lations by Wechsler et al. (2002) . The merger tree model pre- 
dicts slightly lower concentration parameters and a slightly 
broader distribution function. It is important to note that 
both treatments result in quite similar PDFs that are also 
consistent with the results of numerical simulations, despite 
of the competely different assumptions made in each of these 
approaches. 

As mentioned earlier, the uncertainty in the correct 
form of the mass function is quantified by the perturbing 
function G (see equation @). Figure [2] (left upper panel) 
shows how the concentration parameter changes when the 
merger tree is computed with and without the perturbing 
function G. As expected, the concentration parameter tends 
to be larger in the former case, reflecting the earlier forma- 
tion time of halos in the MS as compared with the extended 
Press-Schechter formalism (Wechsler et al. 2002). This re- 
sult illustrates the rather strong dependence of the PDF of c 
on the mass function. This dependence has to be accounted 
for when comparing results from observations and numerical 
simulations. 

The initial conditions of the tree are the concentration 
parameters of the earliest halos. As indicated earlier, the 
particular choice of c(M, z) for the earliest halos does not 
appreciably affect the final value of c, as long as this choice is 
reasonable. For example, Figure [2] (right upper panel) shows 
the probability distributions of c with the initial c(M, z) 
taken from the fit in equation ((6|, and a different fit adapted 
from the results of N07: 



It can be seen that there is no significant change in the 
distribution function. The influence of these initial condition 
on the results is further discussed at the end of section [3731 

The structure of the tree, and hence the calculation of 
c, depends somewhat on the chosen ratio for major mergers, 
q — M>/M<- This dependence is shown in Figure [2] (left 
lower panel); the choice of q is guided by several physical 
considerations. On the one hand, it should not be too small, 
because this would take into account only nearly equal-mass 
mergers. Hydrodynamical simulations (Wik et al. 2008, Mc- 
Carthy et al. 2007) show that mergers with mass ratios as 
high as q = 10 still lead to strong disruption of equilib- 
rium in the inner cluster region, and would thus need to be 
treated as major merger events in our approach. The dy- 
namical impact of taking higher values of q has not been 
explicitly explored in hydrodynamical simulations. Accord- 
ingly, we selected this value to be the highest value of q 
above which the mergers are approximated as continuous 
mass accretion. 

The value of k, which determines the separation at 
which two halos become bound, also influences the results 
quite appreciably. Figure [2] (right lower panel) shows that 
deviations from the fiducial value of k = 5 can shift the dis- 
tribution of c due to changes in cluster initial energies. As 
discussed earlier, the value k = 5 was chosen so as to pro- 
duce realistic distances between clusters that are about to 
merge, and is consistent with estimates of relative velocities 
of merging clusters (Lee and Komatsu 2010). 
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Figure 3. c — M relation at z = (circles), 0.5 (stars) & 1 (trian- 
gles). The thick line shows the results of N07 with the dispersion 
in the logarithm of c indicated by the thin lines. The mass de- 
pendence is consistent with the results of the MS. 

3.3 Scaling relations of the concentration 
parameter 

The expectation values of the distribution functions from 
the previous sections provide the concentration parameter 
averaged over formation histories. It is obviously important 
to follow the redshift evolution of c and its distribution with 
the final mass. Figure [3] shows the c— M relation for several 
redshifts of observation. The results can be well described 
by the scaling relation c oc M~ a , with strong redshift depen- 
dence of a, ranging from a = 0.075 for z — to a = 0.054 
for z = 1, so that the dependence of c on mass is weaker 
for higher redshifts. This likely represents the fact that c de- 
pends on mass through the formation redshift, and the dif- 
ference in formation redshifts for different masses observed 
at z — is larger than for different masses observed at a 
higher redshift. This flattening of the mass-concentration re- 
lation at high redshifts is also seen in numerical simulations 
(Duffy et al. 2008, Gao et al. 2008), although our predictions 
for a are slightly lower at low redshift and slightly higher at 
high redshift than those of Duffy et al. Figure also shows 
the results of N07 for halos at z = extracted from the MS 
(thick line) along with the 1 — a distribution widths (thin 
lines). These results of the MS are consistent with the pre- 
dictions of the merger-tree model, although there seems to 
be a systematic offset between the respective results from 
these two very different studies. 

The dependence on the cluster observation redshift, 
which is often taken to be c ~ (1 + Zobs)^ 1 with 7 = 1, 
is also found to be much weaker and mass-dependent, rang- 
ing from 7 = 0.38 for M = 1O 13 /i _1 M to 7 = 0.24 for 
M — lO 15 /i _1 Af0. This results in slower redshift evolution 
than found by Duffy et al, but is more consistent with the 
findings of Gao et al. for massive halos extracted from the 
MS, especially for masses around M ~ 10 14 AfQ for which 
our result 7 = 0.31 coincides with the evolution seen by Gao 
et al. (note, however, that these authors use the Einasto pro- 
file to describe DM halos). 

In general, the scaling relations deduced from numeri- 
cal simulations are effectively weighted by the mass function, 
and, since the latter has a sharp cutoff at about the typi- 



Figure 4. IC gas temperature PDFs predicted by the merger-tree 
model and best fits to a log-normal distribution for 10 14 /i _1 Mq 
(circles, solid line), and for lO 15 h _1 M (asterisks, dashed line). 
These distributions exhibit a characteristic sharp cutoff at low 
temperatures and a long high-temperature tail, as expected. 

cal mass of a galaxy cluster, mainly reflect the structure of 
smaller, galaxy-sized halos at low redshifts. Extrapolations 
of the results of such simulations cannot faithfully describe 
the structure of massive halos at high redshifts, as pointed 
out by Gao et al. Although we use the results of Bullock 
et al. as the initial conditions for the merger tree - equa- 
tion (J6]) - this choice is justified because our final results are 
not sensitive to the exact form of these initial conditions. 
In addition, the initial halos in the merger tree have smaller 
masses, in the range explored by Bullock et al. 

Full investigation of the c — M relations and their red- 
shift evolution neccesitates the use of numerical simulations 
targeted at massive, cluster-sized halos. We plan to continue 
our study in this direction using the hydrodynamical AMR 
code Enzo. 

3.4 PDF of IC gas temperature 

Since the temperature of IC gas is used as a mass proxy 
in cluster surveys, its PDF is of great observational im- 
portance. We have computed this distribution as outlined 
above. Figure [4] shows the PDFs of the emission-weighted 
temperature for cluster masses M = 10 15 M@ and M = 
10 14 /i _1 M at z = 0. 

As expected, the PDF exhibits a long high-temperature 
tail which corresponds to those clusters that were formed 
atypically early. At low temperatures, on the other hand, 
there is a sharp cutoff that corresponds to clusters that 
formed close to their observation redshift. A log-normal dis- 
tribution provides a good approximation to the temperature 
PDF below M ~ 2 ■ W 15 h~ 1 M Q , as can be seen in Fig- 
ure [4] The width of the distribution is o"io gl0 c = 0.07 for 
1O 15 /i _1 M and cri ogl0 c = 0.08 for 10 14 /i _1 M Q . 

3.5 Temperature scaling relations 

Scaling relations of the gas temperature with cluster red- 
shift, mass, and X-ray luminosity are commonly used in sta- 
tistical analyses of the cluster population and in the use of 
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Figure 5. Mass-temperature scaling relation at z = for the 
merger-tree model (stars) and from observations (circles and tri- 
angles). 



clusters as cosmological probes. Most useful is the T — M 
relation which can be determined from the probability distri- 
bution functions. In Figure[5]we show the emission-weighted 
temperature versus mass for clusters at z = 0. The tempera- 
ture was calculated using equation (|15[) . Blue stars represent 
expectation values of the PDFs, with errorbars indicating 
the distribution variance. The red circles are measurements 
of a sample of clusters from Arnaud, Pointecouteau & Pratt 
(2005), and the black triangles are measurements of another 
sample by Kotov & Vikhlinin (2005), where redshift correct- 
ing factors have been included for both samples. The merger 
tree results are best-fit with the relation T oc M 0,6 , which is 
very close to the theoretical relation obtained for an isother- 
mal sphere, T oc M 2 ^ . It can be seen that the results and 
the expected scatter are consistent with observations. Note 
though the different definitions of mass (M200 in Arnaud et 
al., M500 in Kotov & Vikhlinin) and temperature (spectral 
temperature in both Arnaud et al. and Kotov & Vikhlinin). 

The observational results suggest that the variance of 
the temperature PDF can be seen to represent the amount 
of scatter that is expected in observed clusters due to their 
different formation history. Note that the error in the mea- 
sured temperatures is small compared to the scatter, which 
is slightly larger than the predicted intrinsic scatter, as ex- 
pected, since it has additional contributions. For example, 
not all clusters are fully relaxed and spherical, etc. We find 
that the temperature scales as a power-law in mass at all 
redshifts, T oc M^; however £ varies somewhat with red- 
shift, from 0.6 for z = to 0.63 for 2 = 2, which results 
in slower evolution compared to the simple scaling £ = 2/3. 
We note that the minimal possible temperature of a given 
mass - the low-temperature endpoint of our PDF - scales as 
T min oc M°- 65 ~ - 66 in our model for all redshifts, in much 
better agreement with the standard value. Indeed, the stan- 
dard treatment assumes that the halo is observed immedi- 
ately after it had formed, which is precisely the situation 
described by the low-temperature end of the PDF. The ex- 
pectation value, however, is affected by the width of the 
PDF, which also depends on mass. 

The predicted redshift dependence of T is another key 
relation whose knowledge is important as it reflects on clus- 



Figure 6. X-ray luminosity vs. emission-weighted temperature 
(circles with error bars), where the error bars correspond to the 
disrtibution width in the logarithm of the luminosity and tem- 
perature, respectively. A constant gas mass fraction f g =0.1 was 
assumed. Also shown are measurements of the bolometric lumi- 
nosity vs. spectroscopic temperature from Pratt et al. (crosses). 

ter evolution, and its approximate analytic form is needed 
in comparisons with results of cluster X-ray and SZ surveys. 
The basic redshift scaling of the temperature is contained in 
the relation T oc [E 2 (z) Av (z)] x , where A varies somewhat 
with mass, from A = 0.2 for M = 1O 13 /i _1 M to A = 0.26 
for M = 10 15 /i _1 M©. Thus, T(z) is less steep than in the 
standard relation JTJ, where A = 1/3. In addition, the slope 
of this scaling relation differs with mass, hinting that the 
temperature might not be a separable function in terms of 
mass and redshift. The dependence of these results on the 
model parameters is discussed in Section \3. 81 

The luminosity-temperature relation is an important 
probe of the IC gas. In the framework of the presented ap- 
proach it can be used to test the validity of the simple poly- 
tropic model. Figure [6] shows the luminosity-temperature 
relation obtained from the merger-tree model, as well as X- 
ray measurements of a sample of clusters by Pratt et al. 
(2009). There is reasonable agreement with the data in the 
high-temperature end, with the distribution width approx- 
imately corresponding to the scatter in the measured val- 
ues, but the model clearly overpredicts the luminosity of 
low-temperature clusters. One reason for this could be non- 
constant gas mass fraction which, as hinted by observations, 
is lower in low-mass systems. The dependence of the gas 
mass fraction on mass and redshift could also be related to 
additional physical processes in the IC gas, such as radiative 
cooling and feedback from supernovae and AGN. 

3.6 Integrated Comptonization parameter 

Having determined the IC gas temperature and density pro- 
files (as outlined above), we can now compute another key 
observable - the integrated Comptonization parameter. To 
do so, we also need to specify the gas mass fraction, which 
is taken to be f g = 0.1 for all halos. Figure [7] shows the 
PDF of Y; it exhibits the same general features as the tem- 
perature distribution, a sharp cutoff at low Y , and a long 
exponential tail at high values, largely due to clusters that 
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Figure 7. PDF of the integrated Comptonization parameter from 
the merger-tree model for M = 10 15 /i _1 Mq at z = 0.01. 

formed uncharacteristically early. We note that a log-normal 
distribution is a poor fit to the outcome of our model. 

As in the case of IC gas temperature, scaling relations 
of the mean values of Y can be computed. The scaling with 
mass is YcIa(z) 2 oc M s where cIa(z) is the angular diam- 
eter distance and 8 varies with redshift, from 8 = 1.61 for 
2 = 0.01 to 8 = 1.64 for z = 2. This scaling is close to the 
standard result 8 = 5/3. Similarly, the Comptonization pa- 
rameter scales with redshift as YcIa(z) 2 oc [E 2 (z)A.v(z)] e 
with e = 0.26 for M = 1O 14 /i _1 M and e = 0.31 for 
M = 10 16 /i _1 Mq . The evolution with redshift is slower than 
in the standard description where e = 1/3. 

3.7 Temperature number counts 

The PDFs of cluster observables presented above provide a 
theoretical basis for comparisons with results of cluster sur- 
veys. As an example we consider here the predicted temper- 
ature number counts, which is one of the statistical cluster 
functions that can be used to determine cosmological pa- 
rameters. 

The temperature function, that is the cumulative num- 
ber density of clusters above a certain temperature at a given 
redshift (interval) is computed from the following expression 

n(Ti) = f" hi9h B(Ti\M, z)^f±dM, (20) 

where dn(M, z)/dM is the mass function, namely the num- 
ber of halos per unit comoving volume per unit mass. The 
selection function B(Ti\M,z) is usually defined as B — 1 if 
T(M, z) > Ti and 5 = otherwise, where T(M, z) is found 
according to the standard scaling relations (with a sharp 
cutoff) . 

However, in accord with our treatment here, there is no 
one-to-one correspondence between temperature and mass, 
so we need to incorporate the PDF of the temperature in 
the calculation of the number counts by using the following 
selection function 

f-OO 

B(Ti\M,z)=J P(T\M,z)dT. (21) 

The temperature PDF is described by a log-normal dis- 
tribution with expectation value and variance taken from 
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Figure 8. Temperature number counts: the number density of 
clusters above a certain temperature at z = 0.1 

best fits to the results of our model. The temperature func- 
tions calculated with our more realistic temperature PDF 
and that with the standard relation (between temperature 
and mass) are shown in Figure [8] In the standard calcula- 
tion we chose T(M, z) to equal the expectation value of the 
respective PDF. The calculations were performed over the 
mass range lO 13 -lO 16 ft- 1 M using the Sheth-Tormen mass 
function. 

The two calculations coincide for low temperatures, but 
for high temperatures our improved treatment yields appre- 
ciably higher number counts. The reason for this is that our 
more exact treatment takes into account the long tails of the 
distribution functions. Thus, low-mass clusters with mean 
temperatures below To, that do not contribute to n(To) 
when the standard scaling is used, can have a significant 
overall contribution when the temperature PDF is used. As 
discussed earlier, this is due to the non-zero probability that 
the formation redshifts of these clusters, and hence also their 
temperatures, were higher than the mean values. 

As we have mentioned earlier, the log-normal distribu- 
tion is a mediocre fit to the PDFs of high-mass clusters, and 
a better understanding of their shapes is required in order 
to fully assess their impact on temperature number counts. 
The above calculation demonstrates the importance of tak- 
ing temperature PDFs into account in the analysis of cluster 
surveys. 

3.8 Model uncertainties 

In the previous sections we have shown that our method for 
the determination of the PDFs of the various cluster physical 
parameters provides a relatively simpler procedure to im- 
plement than hydrodynamical simulations. The procedure 
involves specifying several free parameters: q - the maximal 
major merger ratio, k - the parameter that determines the 
initial distance between clusters, and the adiabatic index of 
the gas, r. We should also add to this list the parameters 
of the initial c(M, z) chosen for the smallest halos in the 
tree (see equation (J6j) ) . These parameters were found not to 
influence the results considerably when chosen reasonably, 
in accordance with observational results and N-body simu- 
lations; see the discussion at the end of Section \3. 31 
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of M = W 15 h 1 M but for observation redshift of z = 0.5, 
the following values are obtained: < T >= 12. Ot^ keV and 
s/< (T- < T >) 2 > = 1.6 ±0.3 keV. The results for a mass 
of M = 10 14 /t _1 M© at redshift z = are: < T >= 2.5t ;| 
keV and ^/< (T- < f >) 2 > = 0.5±S;J 9 keV. The relative 
uncertainties in all these cases are similar. 

Finally, we can estimate the robustness of our results 
for the evolution of the scaling relations. As an example, we 
have checked how the scaling of the temperature with mass 
(T oc M c ) and redshift (T oc [E 2 (z)A v (z)] x ) changes when 
the model parameters are varied in the range q = 7 — 13, « = 
4 — 6, r = 1.1 — 1.3. It turns out that £ and A change by no 
more than 1% and 6%, respectively, relative to the values 
obtained in Section VS. 5 1 



Figure 9. Concentration parameter PDF from the merger-tree 
model for M = lO 15 /i _1 M0 at z = and model parameters in 
the range q = 7 — 13 and k = 4 — 6. The thick curve shows the 
PDF calculated with the fiducial values. 



Since we are mainly interested in the global properties 
of the cluster, such as the emission weighted temperature 
and the integrated Comptonization parameter, our results 
have a very weak dependence on a particular choice for the 
IC gas profile. We have repeated our calculations using the 
/^-profile for the gas: 



P( r ) = Pgo 



1+1- 



-3/3/2 



(22) 



with /3 = 2/3. The gas core radius is given by r c — r s /rj 
where r s is the DM scale radius, and a typical value is 77 = 2 
(Ricker and Sarazin, 2001). We have integrated the equation 
of hydrostatic equilibrium to obtain the gas temperature 
profile, setting the pressure to zero at infinity. With this 
profile, the temperature PDFs change by just a few percent 
relative to the polytropic model. We repeated the calculation 
using also the /3-profile with a different boundary condition, 
namely setting the gas temperature at the virial radius to 
the temperature of the IGM, typically 10 6 - 10 7 K. This too 
had only a minor impact on the results. 

In order to estimate the robustness of our model we 
compute the errors on the PDFs that result from small de- 
viations from the fiducial values of the main model param- 
eters. Figure [9] shows the PDF calculated with parameters 
in the range q — 7 — 13, k = 4 — 6. 

Other relevant features of the PDF are the expecta- 
tion value and its variance. These quantities, especially for 
the temperature and Comptonization parameter, need to be 
known in the analysis of cluster surveys. We thus need to 
determine the uncertainty in their values due to variations 
of model parameters. Varying the parameters in the range 
q = 7-13, k = 4-6, T = 1.1-1.3 we obtain < T >= 9.9±°{ 9 4 
keV for M = 10 15 fc~ 1 M B at z = 0, while the width of 
the distribution is y/< (T— < T >) 2 > = 1A±° 2 3 keV. This 
likely is a conservative estimate for the range of the temper- 
ature uncertainty. 

The evolution of the variance of the PDF with mass 
and redshift and its uncertainty can also be estimated. For 
instance, if the parameters are varied in the same range q = 
7—13, k = 4 — 6, r = 1.1 — 1.3, and for the same mass 



4 DISCUSSION 

We have presented an expanded merger-tree treatment for 
the evolution of galaxy clusters that supplements the sta- 
tistical description of the dynamical evolution of DM halos 
with basic physical considerations that enable us to describe 
also the properties of IC gas. It should be stressed again that 
our approach is statistical by construction and is not meant 
to provide a prescription for determining the structure of in- 
dividual halos, but rather to serve as a tool for studying the 
properties of a population of clusters. While our treatment 
is essentially adiabatic, we have adopted an observationally- 
deduced value of the polytropic index. By doing so we partly 
compensate for the fact that gas cooling is not explicitly 
taken into account. Additional justification for the validity 
of our approach is the fact that we are interested here only in 
statistical properties of the cluster population, rather than 
in detailed spatial profiles of the gas density and tempera- 
ture in individual (such as cooling-core) clusters. 

We also assumed that the DM mass profile is not af- 
fected by the IC gas. Although this approximation is often 
made in studies of the statistical properties of a popula- 
tion of clusters (e.g. Bode, Ostriker & Vikhlinin 2009), it 
is likely to be inaccurate when radiative cooling is impor- 
tant, or when there is energy exchange between the DM and 
the gas components, for example during mergers. Numerical 
simulations (Duffy et al. 2010) show that there is a devia- 
tion of at most 15% in the concentration parameter of groups 
and clusters when baryonic physics is included, relative to 
the DM only case. The impact of IC gas cooling on the DM 
density profile is often described by adiabatic contraction 
models (e.g. Gnedin et al. 2004). However, the assumption 
made in these models that the baryons initially trace the 
DM distribution is violated during hierarchical build-up of 
halos. Indeed, Duffy et al. (2010) found that results of the 
simulations were not well described by adiabatic contraction 
models beyond 0.1R V i r . A natural extension of our model 
would be to incorporate IC gas in the halos that constitute 
the merger tree and to follow the joint evolution of both 
components. 

We calculated the PDFs of the cluster concentration 
parameter, its IC gas temperature, and integrated Comp- 
tonization parameter for different masses and redshifts of 
observation. Our deduced PDF of the concentration pa- 
rameter is well fit with a log-normal distribution, in ac- 
cord with results from N-body simulations. The tempera- 
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ture PDF for masses below M ~ 2 • 10 15 /i _1 M© can also 
be described with a log-normal distribution. Our deduced 
mass-observable scaling relations are close to the standard 
relations but contain some corrections - notably the evolu- 
tion of IC gas temperature with redshift is slower than in the 
simple model. The results suggest that the gas temperature 
is not a separable function of mass and redshift. We show a 
possible application of our results to the analysis of cluster 
surveys by calculating IC gas temperature number counts, 
taking into account the effect of cluster formation history. 

The probability density functions of the various ob- 
servables can have important effects on the error estima- 
tion in the analysis of cluster X-ray and SZ surveys. As 
shown by Lima and Hu (2005), large uncertainties in the 
observable-mass distributions may substantially degrade the 
constraints on cosmological parameters from cluster surveys. 
The physically-based estimates of the PDFs of the observ- 
ables considered here provide a tangible basis to begin ad- 
dressing this aspect. 

Among the other related applications of the approach 
presented here a particularly timely one is the calculation of 
the SZ power spectrum, which will be mapped by the Planck 
satellite and several ground-based SZ projects. Comparisons 
of results from our merger-tree approach and those from 
simulations and semi-analytical treatments (e.g., see Sadeh, 
Rephaeli & Silk 2007 and references therein) will yield im- 
portant insight that will help gauging the relative merits and 
disadvantages of these very different approaches. 
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